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ABSTRACT 

Gauss-Hermite and Gauss-Laguerre ("shapelet") decompositions of images have become important 
tools in galaxy modeling, particularly for the purpose of extracting ellipticity and morphological 
information from astronomical data. However, the standard shapelet basis functions cannot compactly 
represent galaxies with high ellipticity or large Sersic index, and the resulting underfitting bias has 
been shown to present a serious challenge for weak-lensing methods based on shapelets. We present 
here a new convolution relation and a compound "multiscale" shapelet basis to address these problems, 
and provide a proof-of-concept demonstration using a small sample of nearby galaxies. 
5P Subject headings: galaxies: structure - gravitational lensing: weak - methods: data analysis - tech- 

niques: image processing 
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1. INTRODUCTION 



With the advent of large photometric surveys, efficient, automated techniques for galaxy modeling have become 
an important part of astronomical data analysis. Models based on a linear combination of Gauss-Hermite or Gauss- 
Laguerre functions (also known as Cartesian or polar "shapelets" , respectively) have become particularly important in 
weak-lensing ellipticity measuremen t, and have also been used f or morpho l ogical classification and in building realistic 
mock image data (see, for 



inst ance, Bernstein fc Jarvis| ( 2002); Refregier] ( 2003[ ); Massey & Refregier (2005); Massey 

^ ( |2004| )T TAndrae et al.Tlpjlol [r~^ 

apelet bases consist ot polynomials weighted by a circular Gaussian, they are a good first- 



etal] ( |2QQ4| ); [Kelly fc McKay| 

Because the standard shapelet bases consist ot polynomials weighted by a circular Gaussian, they 
order approximation to typical galaxy and point-spread function (PSF) morphologies. More importantly, a shapelet 
expansion can be convolved analytically, which substitutes a simple linear algebra operation for what can otherwise 
be the most computationally expensive aspect of a modeling algorithm. 

Recent work has highlighted the drawbacks of standard shapelet-based galaxy modeling, however, and demonst rated 



that even high-order shapelet expansions are often poor representations of real galaxies. 'Melchior e t al.| (20101) have 
shown that these deficiencies can introduce serious systematic ellipticity biases in shapelet-based weak lensmg mea- 
surements. In particular, shapelets cannot reproduce the sharp core and broad wings of galaxies with high Sersic 
indices, and become increasingly distorted at high ellipticities. 

We propose here a shapelet-based modeling technique that can much more compactly represent real galaxies, while 
preserving the lossless analytic convolution and other useful properties of the standard shapelet expansion. By com- 
bining multiple low-order shapelet expansions with different scales, our technique can simultaneously represent the 
extended wings and cuspy cores of real galaxies. We also present a new convolution relation for the ellipse-transformed 
shapelet expansion of Nakajima fc Ber nstein (2007) (hereafter NB07), enabling lossless high-ellipticity shapelet mod- 
eling and eliminating one source ot multiplicative shear bias present in shapelet-based weak tensing methods. 

In Section [2J we summarize the limitations of standard shapelet modeling techniques. In Section |3J we derive an 
exact convolution relation for elliptical shapelets, and we discuss the combination of multiple shapelet expansions into 
a single com pound e xpansion in Section [4J We provide a simple demonstration of the conipound shapelet technique 
using the |Frei et al. ( 1996 ) sample of nearby galaxies in Section [5J and conclude in Section [g] 



2. LIMITATIONS OF STANDARD SHAPELET TECHNIQUES 

2.1. Ellipticity 

The standard Cartesian shapelet basis functions are formed from the product of two one-dimensional Gauss-Hermite 
functions: 



$n(x) = ^ni (^1)^712(^2) 



(1) 

(2) 



This is a clearly an expansion about a circular Gaussian, and is naturally poor at representing functions with high 
ellipticity. However, we can create an elliptical expansion by transforming the coordinate grid by the inverse of the 
transform that maps the unit circle to the desired "basis ellipse" ( |NB07 ): 

(3) 
(4) 
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where a, 6, and are the semi-major axis, semi-minor axis, and position angle of the ell ipse, respect ivelyF] This 



transformation is exact, but the standard shapelet convolution formula (see Refregier 2003| ) applies only to shapelet 
expansions with identical ellipticities; this limits analytic convolution to the case where both the galaxy and the PSF 
are approximately circular. 

Instead, most shapelet techniques make us e of the shapelet-space shear operator, w hich can be represented by a 
matrix multiplication on the basis vector (see Refregier||2003 Massey fc Refregier||2005): 



00 
$„(x|e) = V fSel $„,(x) 

— L J n,m 



(5) 



This relation must be truncated at finite m in practice, however, making the shear operation lossy. A simple elliptical 
Gaussian cannot be represented exactly by a finite circular shapelet expansion, and the approximation introduces 
artifacts that become increasingly severe as the ellipticity increases (Figure [l]). Unless the average galaxy morphology 
mimics these distortion patterns (highly unlikely, given that they involve regions with negative fiux), this necessarily 
makes the goodness of fit of shapelet models worse, on average, as ellipticity increases. While this is most noticeable 



^ An arbitrary rotation can also be applied on the right side, as this corresponds to a rotation of the unit circle; [NB07| rotate by —0 to 
make S symmetric, but we prefer to omit the rotation in order to align the shapelet expansion with the axes of the ellipse rather than the 
original coordinate grid. 
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Fig. 1. — Circular shapelet approximations of increasing order n to an elliptical Gaussian with an axis ratio of 1/4. The number of 
basis functions (|(n + l)(n + 2)) used in each approximation is given in parenthesis. Because the number of basis functions increases so 
dramatically at high order, the order of the expansion must be kept low for computationally feasible techniques, but low-order expansions 
produce distorted representations of even simple high-ellipticity morphologies. 



at high ehipticities, it is important eve n at low ehipticities, as these approximation-induced artifacts systematicahy 
bias shapelet-based lensing techniques (Me lchior et al.|[2QTQ ). 

A better solution is to use ehiptical basis functions to model the galaxy. To convolve with the PSF, one applies the 
inverse shapelet-space shear operator to the PSF to transform it into the coordinate system in which the galaxy is 
round (NB07). However, as the ellipticity of the galaxy model increases, the magnitude of the inverse shear transform 
that must be applied to the PSF model also increases, introducing approximation artifacts in the PSF model. 

This tradeoff is advantageous for two reasons. First, one can often afford to use a higher-order expansion for the 
PSF, because the PSF coefficients are not considered free parameters when fitting an individual galaxy; the PSF model 
is determined separately using images of stars in the same field. In addition, when the galaxy radius is larger than 
the PSF, the undistorted galaxy model plays a larger role in determining the ellipticity of the convolved model than 
the now-distorted PSF model. Still, this technique does not eliminate the "shear artifact" ellipticity bias; it merely 
decreases it by substituting a distorted, approximate PSF model for a distorted, approximate galaxy model. 

2.2. Radial Profiles 

Shapelet expansions also have difficulty reproducing realistic galaxy radial profiles. The azimuthally averaged radial 
profile of a galaxy often follows a Sersic law: 

flux (X e~^ (6) 

with a > 1. Disk galaxies typically have a ~ 1 (an exponential profile), while spheroidal galaxies often have a ^ 4 (the 
de Vaucouleur profile) or greater. The shapelet expansion is based on the Gaussian function {a = 1/2), and hence has 
a much softer core and sharper cutoff at large radii than a Sersic profile with a > 1. In theory, the shapelet basis is 
complete, and can absorb these differences in higher-order terms, but in practice a finite shapelet expansion converges 
to a Sersic model with high a extremely slowly, producing a clear "ringing" pattern in the approximation. Figure [2] 
shows a typical case. 

This poses a clear problem for the use of shapelets as a tool for morphological classification: the slope of the radial 
profile, arguably the clearest and most obvious distinction between galaxy types, is a nonlinear function of high-order 
terms in shapelet space. In contrast, it can be easily be estimated with a simple Sersic fit or concentration measurement. 

The inability to accurately reproduce realistic radial profiles also has implications for shear measurement. Ellipticity 
estimators based on model fitting have been shown to produ ce a multiplicative bias when the model is a systematically 
30or f it to the data (Voigt & Bridle||2010 Bernstein 2010). Using mock galaxies with Sersic profiles, Melchior et al. 



poor n t to tne data ( Voigt 6^ i3ridle||zUlU l^ernstem ^UlU). Usmg mock galaxies witn bersic pronles, Melcnior et al. 
([2010'i have shown that this ''undertitting bias" exists lor shapelets even at relatively high orders. Even at low-surlace 
brightness, modeling the wings of galaxies properly can be very important in shear estimation, a s the wings con tain 
the low-spatial-frequency shape information that is corrupted least by convolution with the PSF (Bernstein 2010). 

A shapelet-like expansi on based on general Sersic profiles rather than Gaussians has been proposed as a possible 
solution to this problem (Ngan et al.|f2009), but these "Sersiclets" lack many of the advantages of the shapelet basis. 



including analytic convolution and last numerical evaluation, and thus far have not been used to construct a practical 
shear or morphological estimator. 

3. CONVOLUTION OF ELLIPTICAL SHAPELETS 

In this section, we derive an exact analytic convolution formula for elliptical shapelets. This provides the "missing 
ingredient" for making use of the elliptical shapelet basis of Equation ^ to address the limitations of the standard 
shapelet basis at high ellipticity. The result allows the PSF model and unconvolved galaxy model to be represented 
by elliptical shapelets with different basis ellipses, and determines an optimal basis ellipse for the convolved shapelet 
expansion that makes the convolution exact. 




Fig. 2. — Shapelet approximations to a de Vaucouleur (Sersic a = 4) profile. The radius and fiux units are arbitrary, and the fit can 
be tuned to perform weh at either small or moderate radii, but not both. Regardless of the scale of shapelet expansion, it will always fall 
off much faster at large radii compared to a pure de Vaucouleur profile. This artificial truncation of the radial profile strongly affects fiux 
and size estimates made using the model, and can bias the measured ellipticity as well, even when the truncation occurs at low surface 
brightness. 



Consider a pair of functions and their convolution, each represented by a finite ehiptical shapelet expansion: 



|n|<Ar/ 
/(X)= J2 /n*n(U-lx) 
n 

\n\<Ng 

n 

\n\<Nh 
(/*5)(X)= Yl /^n$n(W-lx) 



(7) 
(8) 
(9) 



The boldface two-dimensional indices n each run over pairs of one- dimensional indices {ni,n2} as in Equation Q, 
but for conciseness we will define the "magnitude" of a vector index as |n| = ni + n2, as this combination will appear 
often. 

The Fourier transforms of the circular shapelet basis functions are proportional to the basis functions themselves 
(Refregier 2003,), and the Fourier-space ellipse transform is simply the transpose of the inverse of the real-space ellipse 



transform, so 



|n|<Ar/ 

m=\U\ J2 /nil"l*n(U^k) 
n 

\n\<Ng 
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|n|<iV^ 

27rf{k)~g{k) = \W\ ^ hJ^\^^{W^k) 



\p\<Nf\q\<Ng 
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(10) 

(11) 

(12) 
(13) 



We can use the orthogonality of the basis functions to isolate the coefficients of the convolved expansion by multiplying 
by ^m(W^k) and integrating 



\p\<Nf \q\<N,^ 



h^ = 27r\U\\V\ ^ Yl ^'^'+'''"'"' /^'k$p(U^k)$q(V^k)$n(W^k) (14) 

\p\<Nf\q\<Ng 

= 27r\U\\V\ ^ ^ ^l^l+l^l-l"l/p,c,n(U,V,W) (15) 

p q 

To compute the integral, it is useful to separate the exponential and polynomial terms 

/p,q,„(U, V, W) = J(fk<^p{V^k) *q(V^k) *„(W^k) (16) 

= L2]^g-ik-(uu-+vv-+ww-)k ^p(u^k) ^q(V^k) Z„(W^k) (17) 

with 

Zp (U^k) = Up, (Uiixi + U21X2) Up, {U12X1 + U22X2) (18) 

Unix) = (v^2"n!)"'/'i/„(x) (19) 

We prefer the "normalized" Hermite polynomials Tin here over the standard Hermite polynomials H^ both because 
they p rovide a more con cise notation and because their recurrence relations (see Appendix |A| are more numerically 
stable dPress et al.||2QQ2[ ). 



Returning to Equation (17), we can simplify the exponential factor by requiring 

WW^ = UU^ + VV^ (20) 

This reduces to the familiar formula for the convolution of elliptical Gaussians at zeroth order; for an elliptical shapelet 
expansion with ellipse transform S, SS-^ is the covariance matrix of the GaussianP] We then make the substitution 
W^k ^ y: 

/p,q,n(U, V, W) = Jd'ke-^^^'^^^ Zp(U^k) Zq(V^k) Z„(W^k) (21) 

= ^ Jd'ye-y^y Zp(U^W-^y) Z,(V^W-^y) Z„(y) (22) 

Each normalized Hermite polynomial can be represented as a linear combination of monomials: 



N 



nn{x) = Y,Mn,mX'^ . (23) 

m 

where M is the lower-triangular matrix of normalized Hermite polynomial coefficients. We can also write a monomial 
as a linear combination of normalized Hermite polynomials using the inverse matrix: 

AT 

x"^ = Y,M-%Hn{x) . (24) 

n 

^ Note that this relation between the basis ehipses does not assume or require that the elhpticity of the functions we are convolving obey 
the Gaussian convolution rule; while the basis ellipse uniquely determines the elhpticity for the zeroth-order term in a shapelet expansion, 
the elhpticity of a higher-order expansion is also a function of the coefficients. 



Along with the binomial theorem, these can be used to factor the ellipse-transform elements out of the polynomials: 



Z„(Tx) = nn^iTiiXi + T12X2) UnATllXl + T22X2) 
N N-mi 
= 11 Y. Mn„^iM^2,mATllXi+Ti2X2r' {T2lXi+T22X2y 



(25) 
(26) 
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'Hj^{xi)'Hjj,X2) 



(28) 
(29) 



Plugging Equation (29) into Equation (22), we have 

^ |^|<iV/ |m|<Ar, 



/p,q,n(U,V,W) = -- ^ Y. Qp,i(v^U^W-^)Qq,„,(V2V^W- 
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1 m 



|(|<Af^|m|<Af<, 

^ ^ Qp,i(V2U^W-^)Qq,^(y2V^W- 

1 m 

fd^y^i(y/V2J $^(y/\/2) $„(y) 

|i|<Af^|m|<Af<, 
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(30) 



(31) 



(32) 



where Bi^rn,n is the one-dimensional triple product integral defined by Refregier & Bacon (2003); it can be computed 
directly with recurrence relations (see Appendix [b|) . 
The complete formula for convolution of elliptical shapelets is thus 
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S,,,„^„„,(\/2,\/2,l) B^,,„^,,„,(\/2,\/2,l 



(33) 



Because Bi^rn,n{V^^ V^, l) is zero for n > / + m, this relation is exact if Nh > Nf -\- Ng and W is chosen as the 
solution to equation ( [2Q| ), allowing an elliptical shapelet galaxy model to be convolved with an elliptical shapelet PSF 
model without approximation. 

4. MODELING WITH COMPOUND SHAPELETS 

In this section we propose a simple solution to the radial profile problem: instead of increasing the order of a shapelet 
expansion in order to fit Sersic-like profiles, we create a compound shapelet expansion that combines multiple low- 
order shapelet expansions with different radii. By adding a low-order expansion with small radius, a compound basis 
can represent a realistically cuspy core. Likewise, additional expansions with large radii should be able to compactly 
represent extended wings without introducing oscillatory artifacts. 

In particular, we will consider compound basis functions of the form 



j<Np k<N4j] 
j k 



(34) 



While the compound basis thus has a single basis ellipse e that defines the transform Se, it is composed of several 
shapelet bases, each with a different relative scale /3[j]. Each compound basis function may be a linear combination 
of several shapelet basis functions, with weights given by the matrices a[j]k^n' 

A compound basis is thus defined by a sequence of weight matrices ak^n and relative scales f3 (note that all the 
weight matrices will have the same number of columns, but may have differing numbers of rows). When fitting a single 
galaxy, the ellipse parameters e (center, complex ellipticity, and radius) will be considered free parameters, along with 
the coefficients of the compound basis functions. 

Most of the techniques we discuss below also apply to a more general compound basis, in which the component 
shapelet expansions can have entirely different ellipse parameters, but such a basis probably has too much flexibility 
for most galaxy modeling applications. The more general expansion may be of use in PSF modeling, however, especially 
if a physical model of PSF spatial variation suggests that the ellipticity of the PSF varies differently at large and small 
radii. 

4.1. Properties of the Compound Basis 

Two notable properties of the standard shapelet basis are its orthonormality and completeness. Both of these are 
in general broken in the construction of a compound basis. 
The orthonormality condition for standard shapelets 



(l)j{x) (l)k{x) dx = 5jk (35) 

is defined for a continuous integral with infinite limits. While this is an important property for deriving analytic 
formulae involving shapelet basis functions, it is of limited practical use in decomposing image data into shapelets, 
because the pixelization and finite limits of image data do not match the idealiz ed orthonormality c ondition; even 



standard shapelet basis functions are not orthonormal when projected to images (Massey et al. 2004). As a result, 
most standard shapelet decomposition techniques use linear least-squares methods that do not require orthonormality 
of the basis functions, and we will be able to use these same methods with little or no modification in decomposing 
image data into compound shapelets. In deriving analytic formulae for the compound basis functions, it is already 
simpler to start with the analogous formulae for standard shapelets and multiply by the weight matrices a[j]k^n' As 
we will discuss below, however, approximate orthonormality may be useful in some situations, and if necessary, the 
Gram-Schmidt procedure can be applied to a set of compound basis functions to produce an orthonormal compound 
basis. 

The completeness of the standard shapelet basis also rarely useful in practice; while an infinite shapelet expansion 
can perfectly represent any function in R^, in practice we must work with finite expansions. As we have noted, a finite- 
size shapelet basis can be poor at representing the particular functions-those that mimic galaxy and PSF profiles - 
we are most interested in. While the finite-size compound basis thus lacks the completeness properties of the infinite 
standard shapelet basis, our essential argument is that a compound basis can be significantly more complete than a 
similarly-sized standard shapelet basis over the domain of functions we are most interested in. 

4.2. Determining the Basis Ellipse 

For a vector of pixel values z and errors a, we can solve for the best-fit ellipse parameters and basis coefficients in a 
standard least-squares sense: 



min^ — Zm -^^n(^m|e)6n 



(36) 



In general, however, the ellipse parameters e will be highly degenerate with those linear combinations of basis coef- 
ficients b that produce shapelet-space translation, shear, and scaling transforms. For a compound basis composed 
of low-order shapelet expansions, these degeneracies are reduced, because the basis functions poorly approximate the 
ellipse transform. However, even first- or second-order shapelet components will generally make the model too flexible, 
a s first-order terms approxim ate a centroid shift, and second-order terms approximate changes in ellipticity and size. 
Bernstein & Jarvis (2002) propose determining the ellipse parameters using a "null test" algorithm, in which the 
parameters ot the basis ellipse are iterated until the resulting expansion has zero centroid, zero ellipticity, and unit 



size in the ellipse-transformed coordinate system. This is equivalent to solving Equation (36) subject to the constraint 

Y.K^,nbn = Q (37) 
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Fig. 3. — Image, shapelet models (upper), and residuals (lower) for NGC 2683, a typical edge-on spiral galaxy in our test sample. Each 
model has 45 free parameters, equivalent to an 8th-order standard circular shapelet fit. This example clearly illustrates the artifacts 
introduced in circular shapelet fits to high-ellipticity galaxy images; in this case, the circular models actually have negative flux in certain 
regions. As expected, the compound and standard shapelet fits perform similarly at moderate radii, but the compound basis is noticeably 
better at representing structure at small and large scales, modeling the core well while reducing or eliminating the "ringing" residuals 
present in the wings of the standard basis fits. The basis labels are defined in Table [l] and the data and fitting procedure are described in 
section [5] 
with 
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(38) 
(39) 
(40) 
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for some radial weight function W. 

We can construct the analogous constraint matrix for a compound shapelet basis simply by replacing ^ with ^ 
above. Rather than applying the constraint when fitting individual galaxies, however, we can instead "circularize" the 
basis itself by modifying the weight matrices a in equation (34). Given a set of input basis functions ^n(x) and its 
constraint matrix K, we first compute the singular value decomposition of K-^: 

UDV^ = K^ . (43) 

The diagonal matrix D will have at most five non-zero elements, corresponding to the first five columns of U. The 
remaining columns of U give the mapping from the original basis to the circularized basis: 
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(44) 
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(45) 



The resulting basis is naturally five elements smaller than the input basis, and its constraint matrix is the zero matrix, 
so Equation (37) is always true. We can then safely apply the least-squares Equation (36), using the circularized basis. 



to fit both the ellipse parameters and basis coefficients simultaneously. 

4.3. Radii and Shapelet Order 

The primary disadvantage of the compound basis technique is its flexibility. The user must select in advance the 
scales p[j] and shapelet orders N^[j] to be used in fitting an ensemble of galaxy images. This will typically require 
some experimentation on a smaller training sample of galaxies. 



It is not immediately clear what metric to optimize in selecting the scales and shapelet orders, even for an ideal 
training set; one does not want to simply ensure the compound basis fits the mean morphology well. A better choice 
would be to optimize the sum of the x^ values for all galaxies in the training sample, treating the individual galaxy 
parameters as nuisance parameters. This is a huge global nonlinear optimization problem, however, requiring each 
galaxy in the sample to be fit for a new ellipticity and set of coefficients for every iteration in the basis parameters. 
Furthermore, this metric imposes a signal-to- noise weighting on the members of the training sample, which may or 
may not be appropriate for different training sets. Most importantly, it is clear that a basis with a large total number 
of basis functions will generally outperform a basis with fewer, and thus part of the challenge is in determining the 
correct number of basis functions to use. There are also computational advantages to using fewer shapelet expansions, 
just as there are advantanges in decreasing the order of the shapelet components, and a useful metric must assign 
weights to each of these opposing priorities. 

One approach is to set the scales and shapelet orders to match a particular analytic function, such as an exponential 
or other fixed-index Sersic profile. However, to do so, one must choose a metric that defines similarity between the 
analytic profile and its shapelet approximation. While a compound basis can approximate a Sersic function over a 
wider range of scales than a standard shapelet basis, it still cannot reproduce a general Sersic profile over the full 
infinite range of the function, and a metric that puts significant weight at large or very small radii typically produces 
a poor approximation over the range of radii where galaxy profiles can realistically be observed. More importantly, 
even if the analytic function to be matched is a good proxy for the mean or most common galaxy profile, tuning the 
basis to fit a single fiducial profile may reduce the flexibility of the basis, making it worse at fitting morphologies that 
differ significantly from the fiducial profile. We have not explored matching to analytic functions enough to rule this 
approach out as an aid in determining the scales and shapelet orders, but the difficulty in choosing a metric and the 
danger of producing an inflexible basis have caused us to focus our attention on heuristic and training-set approaches 
that do not involve a fiducial profile. 

While we do not have a generalized method to determine optimal values for P[j] and N^[j] for a sample of galaxies, 
a closer look at the form of the first few shapelet basis functions provides insight into how to construct a reasonable 
compound basis. If we consider a typical compound basis as one that contains shapelet expansions at small, medium, 
and large radii, using a zeroth-order expansion for all three will restrict the basis to models with complete elliptical 
symmetry-it will be a sum of Gaussians that share the same center and ellipticity, and differ only in size. As we noted 
in the previous section, first-order shapelet terms approximate a shift in centroid, and as a result adding first-order 
terms to the component shapelet bases of a compound basis allows elliptical isophotes with different centers at different 
radii. Likewise, adding second-order terms allows elliptical isophotes that also vary in ellipticity as a function of radius. 

This is a crucial point in reducing the particular underfitting shear bias noted by [Voigt fc Bridle] (i2010i), as it 
allows the compound basis to represent, for instance, bulge + disk galaxies in which ditterent comiponents have 
different ellipticities. Note that this flexibility is not affected by the circularization procedure of the previous section. 
Circularization does not force the basis functions to have constant ellipticity as a function of radius; it only sets (via 
the weight W) what we define to be the mean ellipticity over all radii. 

While first- and second-order shapelet components can represent some structure beyond the shifted and sheared 
elliptical isophotes discussed here, in general non-elliptical structure will be best fit by adding higher-order terms to 
one or more shapelet components rather than adding additional low-order components. We will explore some of these 
tradeoffs, and the choice of radii, in Section [5) but we caution that these questions should generally be revisited for 
different galaxy samples. 

4.4. Dimensionality Reduction 

Because the radii and shapelet orders must be chosen in advance, a compound basis will typically have a fixed size 
(particularly if it is circularized), and unlike the standard shapelet basis, there is no consistent way to increase or 
decrease the size of the basis to match the signal-to-noise ratio (S/N) and resolution of the image data being modeled. 

Given a training sample, one can attempt to mitigate some of these problems by changing the weight matrices a, 
eliminating those linear combinations of basis functions which do not play an important role in fitting the training 
sample. This can be used to produce a "reduced" compound basis, in which the number of basis functions can be 
much smaller than the number of shapelet basis functions used to build it. Quantifying the "importance" of a basis 
function is difficult, however, and may vary between different modeling applications and datasets. 

A standard application of principal component analysis (PGA) is probably a poor choice. PGA generally requires 
independent and identically distributed measuremen t ve ctors, and neither is generally true for the basis coefficient 



vectors b produced by the fitting algorithm of Section 4.2, unless the training sample has such high S/N and resolution 
that measurement errors can be neglected. In fact, basis coefficients will typically be highly correlated; while the 
orthogonality of the standard shapelet basis results in an approximately diagonal covariance matrix for the fitted 
coefficients, compound basis functions with different radii are not orthogonal, and will produce significant off-diagonal 
terms in the covariance matrix. 

In spite of this, we expect some sort of dimensionality reduction to be a crucial part of any computationally feasible 
weak-lensing or morphological analysis technique based on compound shapelets, and we will explore this more fully 
in a future paper. It may also be an important step in using the training sample to determine the optimal scales and 
shapelet orders of the component expansions; while we have discussed the two problems separately here, constructing 
a basis from a training sample involves optimizing /3[j], N^[j]^ and a^ together. 
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TABLE 1 

Shapelet order at five different radii for the nine demonstration basis sets. 
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Fig. 4. — Image, circular shapelet models (upper), and residuals (lower) for NGC 4123, a typical face-on spiral galaxy in our test sample. 
The elliptical compound basis models the core, bar, and ring structure reasonably well, but none of the models capture much of the spiral 
arm structure. This is generally the case for the face-on spirals in our sample. 

5. DEMONSTRATION AND RESULTS 

To demonstrate the improvements that compound and ehiptical shapelets offer relative to more standa rd shapelet 
techni ques, we apply the compound shapelet modeling method described in the previous section to the |Frei et a~ 
(1996) sample of nearby galaxies. This data set has high S/N and a very small PSF relative to the galaxy sizes, and 
spans a wide range of morphologies. Even though these conditions are not typical of most weak-lensing observations, 
this is an ideal data set for testing the suitability of our basis functions in fitting the true, unconvolved morphologies 
of galaxies, even with application to shear measurement. While the distribution of morphologies for the moderate- 
and high-redshift galaxies used in weak-lensing measurements is not identical to the distribution of morphologies in 
our sample, the same morphological types will be present in both, and a nearby galaxy sample affords relative image 
quality that is significantly better than even the best space-based or adaptive optics imaging of more distant galaxies. 

Of the 113 galaxies in the sample, 82 were observed in R and 5j, and the other 31 in ^, r, and i. We limit our 
attention to the Bj and g images (the two filters are very similar). We fit each of the 113 galaxies in the sample with 
both the elliptical and circular forms of each basis described below, but discard the 31 galaxies in which the fit did 
not converge in one or more cases (these failures are mostly due to large gradients in the background present in some 
of the images). 

We consider nine compound basis sets, including a single-scale standard shapelet basis, each with 45 basis functions 
(before circularization). These are summarized in Table fl] For elliptical shapelet fits (e prefix), we circularize the basis 
with weight function W = \ and fit for the coefficients and ellipse parameters simultaneously using the Levenberg- 
Marquardt nonlinear least-squares optimizer. The initial ellipse parameters are determined from the unweighted 
moments of the image. For circular shapelet fits (c prefix), we fix the centroid and size to the results from the elliptical 
fit, and perform a linear least-squares fit using the uncircularized basis. In each case, there are 45 free parameters. 
We use "forward fitting"Rin all cases to handle the PSF, with a simple circular Gaussian PSF model (with FWHM 
as given in the image headers). 



5.1. Results 



^ The PSF-convolved model is compared with the observed image at every iteration. 
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Fig. 5. — Image, shapelet models (upper), and residuals (lower) for NGC 4636, a typical massive elliptical galaxy in our test sample. The 
difference between the compound basis fits and the standard basis fits is striking, particulary in the oscillating residuals and in the core of 
the galaxy, which is not fit at all by the standard basis. Both circular fits show residuals which are clearly circular, rather than matched 
to the ellipticity of the galaxy, and the standard circular shapelet model is quite clearly circular at large radii, even though the galaxy is 
not; the limited basis simply does not have the modeling power to apply the correct shear at both moderate and large radii. 

The results for three representative galaxies are shown in Figures ^[4j and [5] As expected, the edge-on spiral galaxy 
(Figure |3| demonstrates most clearly the improvement from circular to elliptical shapelets, both for the compound 
basis (4_4_4) and the standard basis (__8__). The difference can also be seen in the c__8__ fit to the elliptical galaxy 
(Figure [5|. Because the standard shapelet expansion already has such difficulty representing the de Vaucouleur profile, 
it cannot afford to "spend" coefficients on matching the ellipticity of the galaxy. The compound basis provides better 
fits to both the inner radii and outer radii of all the galaxies shown, and eliminates the "ringing" residuals except at 
very small radii; this is particularly dramatic for NGC 4636. None of the basis sets provides enough flexibility to fit 
the spiral arms of NGC 4123 (Figure [4| well, and this is generally the case for the face-on spirals in the sample. The 
compound basis fits do provide a reasonably good approximation to the bar and ring structure, and all the fits do a 
moderately good job of fitting the overall profile of the spiral galaxies. 

These tendencies are quantified in Figure |6J which shows the goodness-of-fit distribution for the 82 galaxies suc- 
cessfully fit with all basis sets. The elliptical compound basis consistently has the lowest residuals, and even the 
circular compound basis generally outperforms both the elliptical and circular standard shapelet basis. The elliptical 
standard shapelet basis is not uniformly better than the circular standard shapelet basis, but the cases where it is - 
high-ellipticity edge-on spirals, like NGC 2683 - stand out clearly as spikes in the lower plot, in which the galaxy index 
is consistent across all the fits. It is also worth noting that when the compound basis is used, featureless, low-ellipticity 
early-type galaxies like NGC 4636 have some of the lowest residuals in the sample. In contrast, these have moderately 
large residuals compared to the rest of the sample for the standard shapelet basis fits, demonstrating that using a 
compound basis essentially solves the problem of fitting de Vaucouleur profiles with shapelets. 

Figure [7| similarly shows the distribution goodness-of-fit for all elliptical basis sets. The most obvious result is that 
all the compound shapelet bases tested outperform the standard shapelet basis, both in aggregate and (with one slight 
exception) for every galaxy individually. Between the compound basis sets, those with shapelet expansions at both 
P = 0.25 and P = 4.0 tend to perform slightly better, especially compared to basis sets that contain neither, suggesting 
that compound bases are most effective when constructed with shapelet orders at radii that differ by more than a factor 
of two. This provides a partial answer to the question of how to choose the radii and orders of the component shapelet 
bases that is both frustrating and somewhat comforting: while there is no clear choice for an optimal combination 
of radii and orders, any reasonable choice that samples the radial range of significant morphology is still a marked 
improvement over the standard shapelet basis. 

6. CONCLUSIONS AND FUTURE WORK 

Despite the problems discussed in Section [2J shapelets have become one of the most important tools in galaxy 
modeling, particularly in weak lensing. The inability of the standard shapelet basis to represent morphologies with 
high ellipticity or steep profiles puts limits on the accuracy of shapelet-based shear estimators, however, and likely 
degrades the performance of shapelet-based morphological analyses. We have presented here two techniques to address 
these limitations. 

The first is a convolution relation for elliptical shapelets. This eliminates the need to use the lossy shapelet-space 
shear operator, and allows ellipse-parameterized shapelets to be used in place of a circular basis parameterized only 
by centroid and radius. While the elliptical basis requires five nonlinear parameters to be fit instead of three, the two 
additional parameters more than pull their weight in modeling power, particularly for high-ellipticity galaxies. Perhaps 
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Fig. 6. — Rank-ordered goodness of fit for 82 galaxies from the |Frei et al.| |l996| sample of nearby galaxies, for selected elliptical and 
circular shapelet basis sets. In the upper plot, the galaxies are rank-ordered by reduced x^ for each basis separately; each point on the 
X-axis thus corresponds to the nth best-fit galaxy for each basis. In the lower plot, the reduced x^ from the e__8__ fit is used to set the 
galaxy rank index for all bases; a point on the x-axis thus corresponds uniquely to a single galaxy. The goodness of fit for the galaxies 
shown in Figuresjs]^ and [s] are given by the square, circle, and triangle points, respectively. The circular standard basis actually produces 
slightly lower average and worst-case residuals than the elliptical standard basis, but the compound elliptical basis is significantly better 
than either. The high-ellipticity edge-on spirals (such as NGC 2683) that present a particular problem for the circular shapelet expansion 
can be clearly seen in the lower plot as spikes in the circular basis goodness-of-fit. 
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Fig. 7. — Goodness-of-fit comparison for all elliptical basis sets. The rank-ordered indexing and galaxy sample are the same as in Figurelo] 
but the y-axis scale has been changed to show more detail. All of the compound shapelet basis sets outperform the standard shapelet basis 
for all galaxies but one, in which the results are comparable, and the compound bases show considerably better worst-case performance. 
Between compound basis sets, those with components at both the largest and smallest radii (4_4_4 and 206_3) consistently have the lowest 
residuals. 
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most importantly, an exact elliptical shapelet convolution allows shapelet-based shear estimators to be constructed 
that do not suffer from the "shear artifact" bias discussed in Section [2TT1 The results of Section |5] do not demonstrate 
the advantages of the elliptical convolution formula; while we used it in fitting the galaxies, the size o f our te st galaxies 
relative to the PSF makes the difference between our exact convolution and the approximation of NB07 negligible. 
However, this formula can be immediately put to use in existing elliptical shapelet shear-measurement methods, such 
as that of NB07, eliminating one source of bias for galaxies near the resolution limit. 

To address the difficulties of the standard shapelet basis in fitting galaxies with high Sersic indices, we have introduced 
the concept of a "compound" shapelet basis - a basis composed of multiple low-order shapelet expansions with 
different scale radii. Even with simple, ad hoc choices for the radii and orders of the component shapelet expansions, 
the compound basis is significantly better than a single higher-order shapelet expansion at fitting realistic galaxy 
morphologies, particularly in its worst-case performance. By combining low-order shapelet basis functions at multiple 
radii, a compound basis can compactly represent both the sharp cores and extended wings of galaxies with high Sersic 
indices. As we have seen, this can make a significant difference even in spiral galaxies with relatively low Sersic indices, 
and almost completely eliminates the oscillatory artifacts often present in standard shapelet fits. This should reduce 
the underfitting shear bias for shapelet-modeling shear measurement methods, but whether the bias can be reduced to 
acceptabl e levels for future su rveys is a qu estion we reserve for a future paper. Morphological analyses such as those 
of Kelly & McKay (2004) and Andrae et al. (2010) should also benefit from utilizing compound shapelet basis sets, as 
the more compact compound shapelet representation essentially gives classification and clustering algorithms a head 
start in reducing the dimensionality of the problem. 

Unlike the elliptical convolution formula, further work is needed to make full use of the compound shapelet concept 
in a complete shear measurement or morphological analysis method. While we have demonstrated success in fitting 
nearby galaxies with unoptimized basis sets and a simple fitting algorithm, more challenging modeling problems may 
require techniques that tune the parameters and weights of a compound basis using a training sample of galaxies. 
Such a basis of "eigenmorphologies" , built from analytically tractable shapelet building blocks, would be an extremely 
powerful tool, not only for morphological analysis, but also in further quantifying and reducing underfitting biases in 
model-based shear measurement techniques. 



The author thanks Tony Tyson for very helpful comments on earlier drafts of this paper, without which several 
figures would have been nearly incomprehensible. 
This material is based upon work supported under a National Science Foundation Graduate Research Fellowship. 
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APPENDIX 

A. NORMALIZED HERMITE POLYNOMIALS 

The normalized Hermite polynomials are simply the product of the standard Hermite polynomials with the Gauss- 
Hermite normalization factor 



7r2^n!) ^^^ Hn{x) . 



Like the standard Hermite polynomials, these can be efficiently evaluated using recurrence relations 
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These recurrence relations are more numerically stable than those for the standard Hermite polynomials, which suffer 
from round-off error at high order. 

B. TRIPLE PRODUCT INTEGRAL 
In Section [3J we defined the elliptical shapelet convolution formula ( [32| ) in terms of the triple product integral 
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with Q/ = /3 = >/2 and 7 = 1. The recurrence relation given by Refregier fc Bacon| (2003) to compute this integral 
suffers from the same round-off error problems as the standard Hermite polynomial recurrence relations. Using the 
normalized Hermite polynomials, we present here a more stable relation. Let 
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Note that while our B is the same as that defined by [Refregier fc Bacon (2003||, our C is not the same as their L: 
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